General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 




AERONAUTICAL 
AND ASTRONAUTICAL 

ENGINEERING DEPARTMENT 




(NASA-CE-153913) DESIGN OF HIGH LIFT N77-27108 

AIRFOILS WITH A STRATFORD DISTRIBUTION BY 

THE EFPLER HETHOO (Illinois UniT.) 86 p HC 

A05/HF A01 CSCL 01A Unclas 

G3/05 39191 


to NASA STI FACIUIY S'] 
input BRANCH c?/ 




, ‘ ‘ ENGINEERING EXPERIMENT STATION, COLLEGE OF ENGINEERING, UNIVERSITY OF ILLINOIS, URBANA 


I 



Aeronautical and Astronautical Engineering Department 
University o£ Illinois Urbana, Illinois 


Technical Report AAE 75-5 
UILU-Eng 75 0505 

NASA Grant NGR 14-005-144 
Allen I. Ormsbee, Principal Investigator 


DESIGN OF HIGH LIFT AIRFOILS WITH A STRATFORD DISTRIBUTION 
BY THE EPPLER METHOD 

by 

William G. Thomson 


University o£ Illinois 
Urbana, Illinois 
June 1975 


TABLE OF CONTENTS 


CHAPTER PAGE 

I. INTRODUCTION 1 

II . THE EPPLER METHOD. ... ... 3 

III. THE STRATFORD DISTRIBUTION. 20 

IV. USING THE EPPLER METHOD TO DEVELOP A STRATFORD 

AIRFOIL. . . ........... 33 

V. SUGGESTIONS FOR FURTHER STUDY. 48 

APPENDIX 

A. THE EPPLER PROGRAM. ... .... 55 

B. THE STRATFORD PROGRAM ............ 66 

C. THE LOCKHEED PROGRAM................................. 69 

LIST OF REFERENCES. .... .... ................... ....... 73 


I. INTRODUCTION 


High lift airfoils have been the subject of on-going study at the 
University of Illinois [l,2,3,4,5] for several years. The major part 
of the present study has been concentrated on airfoils having a Stratford 
[3,6] pressure distribution, which has zero skin friction in the pres- 
sure recovery area. This pressure recovery represents an approximation 
to maximum pressure recovery without separation. 

These airfoils are designed with non- zero velocity at the trailing 
edge. This non-zero trailing edge velocity is unavoidable, as the zero 
skin friction velocity approaches a zero velocity tangentially, and it 
would require an airfoil with an infinitely long chord to have a stagna- 
tion point at the trailing edge if a Stratford distribution was used. 

Chen [3] determined the optimum relationship between the maximum velocity 
on the upper surface (the "rooftop" velocity) and the trailing edge 
velocity, and this relationship has been used since then in the design 
of high lift airfoils at the University of Illinois. However, this 
relationship does not specify the magnitude of either the rooftop 
Velocity or the trailing edge velocity, but only specifies the ratio 
between the two velocities. Therefore, as the trailing edge velocity 
increases, the lift increases (since the rooftop velocity goes up) , 
and the maximum lift possible is limited only by the fact that, as the 
trailing edge velocity increases, the angle of the trailing edge in- 
creases, distorting the shape of the trailing edge. 

The Eppler [7,8] program is an inverse conformal mapping technique, 
where the x and y coordinates of the airfoil are developed from a given 
velocity distribution. Unfortunately, the velocity distribution is ■ 
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given in terms of the circle plane, and the transformation from the circle 
plane to the airfoil is not known until the problem is Solved, so the 
Stratford distribution cannot be used as a direct input to the Eppler 
program. The problem is further complicated by the fact that the roof- 
top velocity is not known until the problem is solved, so the desired 
velocity distribution (the Stratford distribution) is not known until 
after the shape of the airfoil has been determined by the Eppler pro- 
gram. Therefore, the solution of the problem involves visually compar- 
ing the output of the Eppler program with a Stratford distribution, and 
then guessing the modifications to the input of the Eppler program to 
get the desired output. With experience in determining the changes re- 
quired in the input to the Eppler program to yield the desired changes 
in the output, the number of iterations required to determine the correct 
airfoil decreases to a reasonable amount. 
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II, THE EPPLER METHOD 


The method developed by Eppler 17,8] is an inverse conformal map- 
ping technique that determines the x and y coordinates from a given 
velocity distribution. The two planes involved are shown in figure 1, 
The C plane shows the flow about a circular cylinder, while the z plane 
represents the flow about the airfoil. The velocity in the z plane is 
given in terms of coordinates determined in the ^ plane, z and ? are 
defined as: 

z = X + iy (1) 

s = 5 + m = re^^ C2) 

The flow in the ^ plane is such that the rear stagnation point falls on 
the real axis at ^ = 1 . 

There exists a transformation of the ? plane to the z plane such 

that the z plane represents parallel flow about a closed airfoil at an 

angle of attack a. Since C = 1 represents a stagnation point, the Kutta 

condition requires that this must transform to the trailing edge of the 

dz 

airfoil. As this is to be an infinite parallel flow, z(oo) = «> and 
must be real. The general function that satisfies these requirements is 

oo 

K) . ft Z * J & (3) 

v=o " ■ 

where is real, but not equal to zero. 

The complex potential in the ^ plane can be represented as 

F CO = (J) + i'i' = CCe'^“? ^ In ^ C4) 



Figure 1. The complex mapping planes 
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where P is given by 


r = 47 tC sina 


C5) 


The complex velocity in the z plane is given by 


Ve 


■ i6 _ ^ dF/dg 

dz dz/dC 


C6) 


The inverse of this will be used, or 


dz dz/dg 
dF ■ dF/dg 


(7) 


In order to prevent an undesirable root, this can be written as 


, dz ' dz , dF 
ap = dc ■ d? 


( 8 ) 


The velocity vector in the z plane can be introduced as V 


\r 

Ve . 


Then 


dF ,, -i0 
T- = Ve 
dz 


C9) 


Therefore, 


In ^ = - In V + i9 
dF 


( 10 ) 


d z 

The real part of In is then -In V. Outside the boundary of the unit 
dz' 

circle, Iri is regular, and can be calculated when the real part is 

known on the boundary. Since F(g) is also known (eq.(4)), equation (8) 
dz 


can be solved for z (g) then can be found by integration. 


2 (O 


must be of the form of equation (3) to match the boundary conditions . 

dz 
dF 

in a satisfactory 


Therefore, the problem is to find an equation for In which results 
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From equations (4) and (5) 


dF ia ^ r 1 ^ ^"2ia, 

dc" C - 1) ( J * a ) 


( 11 ) 


In light of the singularities involved at the stagnation points. 


dz 


In ^ can be represented by 


dz ^ r_ici y 1 , _-2ia^j y + ib ) e"’" (12) 


In-^ = T^ln C - In [e^^ ( ^ + e 

ur (, 


m=0 


Using this equation and equations (8) and (11), 

GO 

in (1 - h + J: Ca„ * ib„) C‘” 


a 

This results in 


r' m in 

^ m=0 


dz 1. 

di; "" ® 


y (a + ib ) C 
m=0 


-m 


(13) 


(14) 


Expanding this yields 


iz ,, 1, % ^‘‘o (^2 * 

jj= Cl - 5) a e e 


(15) 


d z 

z(5) must be of the form of equation (3), so must be of the form 


dz o n V .0 >.-V-l 


V=1 


(16) 


If we let 5= + ib^, equation (15) can be further expanded to 


Cl - h Cl VA + ^ ♦ ...) (1 * *1* t *1* 

= 21,2 31^ 

C17) 


^ ^2 

+ .,4 ( 1 +—-^+ 


■" y~ 
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Comparing equations (16) and (17) yields 


Pi = 1 * \ " fr * iV " 


. = e 


(18) 


d z 

Since 0^ is real, is also real. At infinity, we want = 1. indicat- 
ing that the flow at infinity in the circle plane is in the same direction 

dz 

as the flow at infinity in the airfoil plane. From equation (16), 
at infinity. Therefore, 

A 


0^ = e " = 1 


(19) 


Because of this, A =0. 

o 

Comparing the terra yields 

^1 ^1 ^1 
_ =0 
c c • 


( 20 ) 


So a^^ = 1 , = 0. 

At infinity, we can arbitrarily set the velocity equal to unity. 

From equation (12), in the limit as ~ = Ce"^*^. Therefore, C=1 

and In C=0. 

The problem that remains is to define the a and b that have not yet 

ni m ^ 

been defined such that, along the surface of the airfoil, the velocity 
assumes the prescribed values. 

Along the surface of the circle plane, C - Using this,j 

equation (12) results in 


In = -In (e^"^ + e^^°^ ) ] + I (a + ib ) 

dF 1. ^ j j L m ra-^ 

m=u 


( 21 ) 


If we use the substitutions 


P((j>) - I (a cos m^ + b^^ sin m(|)) 
. . ro=0 
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QC<jj) = I sin nuf) + cos mtj)) 


(23) 


m=0 


equation (21) can be written as 




+ iQ((j)) 


(24) 


This can be written as 


-In V(^) + i6((j)) = -In 2|cos i(^/2-a) ]+ ?(())) + + Q('j>) 

- {{tt}}} 


(25) 


where {{tt}} is given as 


° (Ql^XT^ + 2ot) 
^ TT (iT+2a<c|)<2Tr) 


(26) 


The {{it}} term is necessary due to the shift in the direction of the 
velocity at the stagnation point. 

The real part of equation (25) can be rearranged into the form 


P((j)) = In 2 1 cos C‘|- - a) j -ln V((j)) 


(27) 


Through harmonic analysis, the a^^'s and b^’s can be determined 
from equation (27). However, we must have a^ = 0, a^ = 1, and b^ = 0, 
due to equations (19) and (20) . Therefore, 

2tt 


P(^) dc|) = 0 


(28) 


2it 


p (4>) COS^ d(|) = IT 


( 29 ) 


2ir 


P((|>) sin^ d^ = 0 
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The velocity distribution we specify must meet these requirements. 

By integrating equation (14) , the transformation 

00 




(1-1/0 e 


I Ca. 

m=0 


+ ib ) 
m m-' 


.-m 


dC 


(31) 


id) 

can be derived. This yields the flow for the entire z plane. If t,= e 
is entered into this transformation, the resulting z = x + iy will yield 
the profile of the airfoil . The results are 


x(4>) = I -4 sin I I cos (|- - a) | cos (^ + Q((j))) dcj) (32) 

■ y(0 = I -4 sin | | cos (| - a) | sin (-| + Q(tJ))) # (33) 


The only quantity remaining to be defined, then, is Q((|)) . However, 
Q((^) is a conjugate harmonic function of P(0> and can be derived from 
the formula 


27r 


QC4>) ■= 


_i_ 

2tt 


cr 


P(o) cot 




da 


(34) 


Given a velocity distribution that yields a P(cj)) such that equa- 
tions (28) through (30) are satisfied, and an angle of attack, the x and 
y coordinates of the desired airfoil can be generated using equations 
(32) and (33). The angle of attack, a, need not be held constant, but 
can be a function of (j). Thus, the upper surface can be designed at a 
different (higher for reasonable airfoils) angle of attack than the low- 
er surface, or even different portions of the upper (or lower) surface 
can be designed at different angles of attack. 

The profile of the airfoil is determined by a and b . Therefore, 
fo.r a fixed p rof ile, and b^^ are fi xed. Altering the angle of attack 
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will not alter the airfoil profile, and, therefore, will not alter a^^ 
or b^. This means that P((|)) is independent of the angle of attack. 
Equation (27) can, then, be written in the form of 


P (cf>) = In 2 1 cos ('I' - a ) j - In V ((f)) 


where V (^) is the specified velocity at the point on the airfoil cor- 

'h 

responding to tj), and a (<J)) is the corresponding angle of attack. At 
any angle of attack, a, then, the velocity can be given as 


v*(« 


The circle plane can be divided into segments, as in figure (2), 

where . . . • (jJj • <|>j indicates the stagnation point. 

L S' L 

The angle of attack specification takes the form 


a ((()) = = constant for 


and the velocity takes the form 


V C(|)) = V.W. (<|)) 


where is a constant for (b^ and W((|)) is given as 


cos(b - cosij), , - 1 -p 
^(0 


“1 -Pr CO 

} 1 - 0 . 36 {{— 


cos6 - cosi 2 -1 H 


on the upper surface i On the lower surface, the velocity distri- 
bution is similar, but different values of IC,, p, (j) , and are used, 

n w S 

indicated by 1^, p, and respectively. The terms in the double 
brackets of equation (38) are defined by 

tf(^) (f(^) > 0) 

{{f(<j>)}} = C ■ C40) 

{ 0 (f(« < 0) 
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Equation (39) can be considered to be of the form 

K 

W(<j)) = W^((J)) W^((j)) “ (41) 

or, on the lower surface 

- - 

W(c|)) - W^(4>) (42) 

The (or W^) term produces the major pressure recovery. This is 

in the form specified by Wortmann {9], in which the shape parameter is 

held constant, which delays separation. The (or W^) term develops 

the cusp distribution. It is generally applied to the last 3-5% of the 

airfoil length. Outside of the range of the specified region, for 

W or for W (or or (J) on the lower surface) , the value of 

w s s w s 

W or W is as given by figure (3). 
w . s 

P((jO must be continuous over the airfoil, so P(^^") = P(4>j^^)» 

* * 

and P(0) = P(2 tt). Substituting our values for V and a into equation 
(27) yields 

K 

P((j)) = In 2|cos (|- a((j)))| -In “ (43) 

Since, at the trailing edge, P(0) - P(27r), 

K 

ln2lcosa. |-ln [V,W^(0)W = ln2 j cosaj | -In [Vj. iT (27r) 

a. a 

- 

W^(2ir}] (44) 

At all other segment boundaries, which are generally outside of 
the cusp region, = 1, P(<Ji^ ) = PC<J>j^^) » and therefore 

ln2|cos(-^ - olj) -ln[V.W^(y] = ln2|cosf l 
: --In ' ^ 
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On the lower surface, is replaced by in equation (45) . 

If, in any segment, (p = tt + 2a^, P((p) becomes "infinite. This 
is most likely to occur in the segments to either side of the stagna- 
tion point, + 1 » where <f).. indicates the stagnation point. 

■^L* -L L 

In order to prevent this, we require Cy 


a > ot 

h 


(463 


and 


7T + 2a^ > (pj > TT + 2a^ . 

L L 


(47) 


If all the (|)7s and a^'s along with y, U» and <p^ are 

given in the problem specification, this leaves only and left to 

be determined. However, with equations (28) through (30) there are three 

conditions that must be met. Therefore, we need to free one of the 

quantities listed above. The quantity best suited for this is tf>_ , 

. ‘ . L 

the location of the stagnation point. 

If we use equation (43) to define P(<f>), equation (29) becomes 

' 2tt 


0 


[ln|cos(-|- - a((J>) j -In V((|)) -In W^(<j))-K^lnW^ (({)) + lu2] cos(f)d^=7r 

^■■■''.'■(48) ■ 


This can be expanded to 
I 


>i 


[ [In] cos (y - a. ) I -InV . -InW -K„lnW +ln2] cos(|)d(i) 

• ^ 1 J * ^ ^ W n S 

A X ... 


h-1 

I ({) . 

J 

, h-i 


[In I cos - a^) |-lnV^-lnW^-K^lnWg+ln2] cos(J)d(|)=7r 


( 49 ) 


csaasasa»6taaaBf3?sais:a 


■--m 


The integral 


cos(j) In j cos ^ - ot^|d({) can be evaluated 


as 


coscfilnjcos ^ - oi^|d(f) = (sincji + sin2a^) In | cos '2 “ *^3^ 
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1 C 50 ) 

+ -j Ctl^cos 2di^ - sin<j)) + constant 


If we denote W . as 
cx 


W . = - 
cx 


(j). 


•'i-l 


cos(|) In ((j)) dc|> 


( 51 ) 


and introduce the notation 


4 >i 

In (i.j) = Injcos (-^ - a.) 


(52) 


equation ( 49 ) can be written as 

I 


^cl ^ ''^cl ^ ^^i (i-l>i)) 

a i=l 


1 1 V 

+ -9 (ff)^ “ 1) cos 2a. + 95- (sincf). - sincj). y) + sin(|). tln(i,i) 

^ 1 j-“l X Z 1-i X 


■InV. ] - sin<j) [ln(i-l,x) -In V.]} - 

X X 

r2TT ^ 


cos(p lnW^d(j) 


cos(|) InW^'(^) dcf) = 7T 


(53) 




w 


Due to equation ( 45 ) , the next to the last term, with i=n, and the last 

term, with i=n+l (n=l, I - 1 ) , cancel each other out. Out of these terms, 

«• , 

only the last with i=l and the next to the last with i=I remain. However, 

a ■ 

since (J)^=0 and cf)j. = 2'ir, sin$^= sintj)j = 0, these terms also drop out. 
a ^a 

The third term of the summation also drops out. Therefore, equation ( 53 ) 
can be written as 
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‘'^cl * ^cl -ln(i-l,i)) 

a i=l 


+ -J C(f>^ ~ (f)^_j) cos 2a^} - 
2rr 


(p 

r^w 


cosi In W d(b 

yj r 


(p 


cos(p In (<p) d(p = TT 


(54) 


w 


We now introduce J such that 
c 


a 

2a>((lnCi4)-ln (i-l,i)) +4 (<#>. -(p:, J cos 2aJ 

^ l l-J. X 


w 


2ir 


0-' 


cos(P In W dd> - 
w ■ 


<P^ 


cos(p In W^ d(p 


w 


If we define a = ot- . .= 0, We can alter the indices to get 

O 1 + -L 

I ^ cP 


y 

- TT = > 


{sin 2a. - cos 2a. ,)} - 
1=1 1 1+1 


cos(j) In d(|) 


(55) 


2nV 


<P. 


Cosfp In d(p 


w 


Now equation (54) can be written as 


(56) 


K W , + K„ WpT + J - TT = 0 
H cl H Cl c 
a 


(57) 


In a similar manner, if we define W . as 

SI 


W . = 

SI 


sin(^ In W C^) d(J> 




C58) 


and we define J as 
s 


= I -d+cos 2a^.) In (i,i) + (1+cos 2a^^j^)ln Ci4 


i=0 


+ 2 ^ I (cos 2a^ - cos 2aj.^j) 


<P^ 


w 


2TT 

f 


C0S(J> In dc|) - J cos^ In d<j> 




w (59) 


o 
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Equation (30) can now be written as 


K^Wgj.KaWsl *^3=“ 

a 


( 60 ) 


From equation (45) , 


In = In (i,i) -In (i,i+l) + In 


(61) 


Substituting tliis value o£ into equation (44) yields 


In Icosaj I “In (i,i) + In (i,i+l) + In In W^^(O) 


= In [cosa^ l“ln In (2 tt) 


(62) 


From equation (45) , V can be determined as a function of V™, and so 
forth until is reached. In this manner, and Vj are eliminated 


from (44) yielding 


00 


-K^ In W^(0) + Tpj In (271) + I {-In (i,i) + In (i,i+l)} = 0 

(63) 

Defining as being the summation term of this equation, equation (63) 
can be written as 


■"h 


In W^(0) + Ky In (2 tt) + = 0 


(64) 


We now have three equations (equations (57) , (60) , and (64) ) for 


the three unknowns (Kj^, Kj^ and tj)^ ) . Eliminating and yields 

- L . 


CJc - tt) Dj - . O 


(65) 


where 


“l = «sl + *«SI 1" \ ™ 

3. 


( 66 ) 
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°2 = ^cl * ^CI 

a 


D„ = W , W ^ - W ^ W^. 

3 cl si cl„ si 

a a 


(67) 

( 68 ) 


Equation (65) is a transcendental equation for c|)^ . Once the value of 

■‘•L 

(j)j is determined, and Kj^ can be determined from equations (57) and 
(60) . 

With, equation (45), we now have I - 1 equations for the I values 
of V^. The last equation comes from the previously unused equation (29), 
which guarantees uniform flow at infinity. 


2tt 

P((j)) dd) 


0 


I { 

i=l 


^i 

r 


[In I cos (|- - a^) 




V. 

■In 

1 


i-1 


In - In Wg(d)) ] d<t>} + 27r(ln 2 - In = 0 

(69) 

Now that all of the Kj^, and are known, P(<j)) can be calculated 

from equation (43) , Q((j>) can be calculated from equation (34) and x((j)) 
y({{j) can be calculated from equations (32) and 03) . 

For practical numerical calculations, the circle plane is divided 
into 2N equal parts, with the positions given by 




vir 

N 


(V = 0, 1, 2, .. ., 2N-1) 


(70) 


Next, the d>. 's (except 4 ) a. 's, K, K, y, and y are chosen. The values 

■ H > 

of W^J, Wgj can then be calculated (equations (51) and (58)). 

. ; a.' ^a. 

Using equations (66) through (68) , D 2 , and can be calculated, 

and the transcendental equation can be established. Once is dete.rmined 

■ 
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and can be calculated from equations (57) and (60). P((j)) and 

Q((j)) can be determined at each point on the circle. Then x((f>) and y(4)) 

are determined, so 2N points are determined on the airfoil. These 

points are equally spaced on the circle plane, but they are not equally 

spaced in the airfoil plane. *■ 

The resulting coordinates will yield an airfoil oriented at 

its zero lift line. However, the angle between the zero lift line 

and the chord line (3) can be dv^terminsd and subtracted from a to 

yield the angle of attack with respect to the chord line. 

Since the values of and K , are determined by the closure re- 

n n 

quirements, there is no direct control of these values in the input 

specifications. and determine the trailing edge angle. In 

order to maintain some control over this trailing edge angle, a 

desired value of can be specified, and by iteration, 

varying either on the upper or lower surface (or both surfaces), 

or K or K, the desired K can be attained. 

s 
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ni. THE STRATFORD DISTRIBUTION 


The problem of high lift airfoils with a Stratford distribution 
was investigated by Ormsbee and Chen [2,3] in 1971. The lift of an 
airfoil is generally limited by boundary layer separation on the upper 
surface, where the fluid is subjected to an adverse pressure gradient. 

The usual velocity distribution on an airfoil generating lift is one 
in which the flow accelerates as it goes around the leading edge, reaches 
a maximum, and then decelerates as it approaches the trailing edge. The 
problem of separation occurs during this deceleration of the flow. 

In order to achieve the maximum lift, we want to accelerate the 
flow on the upper surface to a maximum, hold this maximum velocity 
for as long as possible, and then decelerate the flow as quickly as 
possible without separation. This is similar to the problem investigated 
by Stratford in 1957 [6]. Stratford investigated flow over a flat plate, 
decelerating the flow by varying the shape of the v/all facing the test 
surfaced On an airfoil, the velocity gradient is developed by varying 
the shape of the airfoil, but the results should be the same no matter 
how the velocity gradient is obtained. 

In developing the Stratford distribution, incompressible flow at 
a large Reynolds number is assumed. The lift, then, is given by the 
Kutta-Jowkowski theorem. 


L = p U^ r i (71) 

Therefore, assuming a fixed density and free stream velocity, the only 
way to increase the lift is to increase the circulation. The Circulation 


r, is defined by 
s , 


r = 


V(s3 ds 


(72) 
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where s is the distance from the trailing edge to the point in question 
measured clockwise along the surface of the airfoil, and s^ is the total 
length along the surface of the profile of the airfoil. V(s) is the local 
velocity on the surface of the airfoil. If we use the subscript L to 
indicate the stagnation point, the integral for T can be broken into 


two separate integrals 


r = 


V(s)ds + 


V (s)ds 


(73) 


Along the lower surface, from s = 0 to s = the velocity is always 

^ fS. 


in the direction of negative s. Therefore, 


V(s)ds is negative. 


. 0 . 


and the largest value possible is zero. This indicates the velocity 
along the lower surface is always zero, or stagnation occurs along 
the entire lower surface. Assuming zero velocity along the lower sur- 


face, the circulation can be given as 




V(s)ds 


(74) 


For convenience, the origin of s may be shifted from the trailing 
edge to the front stagnation point, and the distance along the upper 
surface may be referred to as s^, such that 


s = s„ - s, 
u T L 


(75) 


In this case, equation 
s 


may be written as 


u 


r = 


V (s) ds 


(76) 
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The lift can be normalized with respect to the free stream dynamic 
head and the length of the upper surface. 


f = “ 

L 1 „ 2 

TT D U S 
2 ^ u 


and the velocity can be normalized with respect to the free stream velocity 


qCsJ = f-) 


Equation (71) can then be written as 


s 


q(s)ds 


As mentioned earlier, the velocity q(s) starts out as zero at 
the stagnation point, increases rapidly to a maximum velocity, then, 
starting at a point s = s^, decreases as rapidly as possible without 
separation down to the trailing edge velocity. Boundary layer separa- 
tion is imminent when the local skin friction goes to zero. Therefore, 
in order to have a maximum deceleration, we want to have a velocity 
gradient such that the skin friction is zero from s = s^ to s - s^. 

Stratford derived an expression for flow over a flat plate with 

zero skin friction based on C , the pressure coefficient based on 

PQ ■ 

the pressure before the pressure rise begins. 


: ^ - Pt 

>0 ‘ Po - Pt 




where is the total pressure, p^ is the static pressure before the 
pressure rise, p is the local pressure, and is the initial velocity, 

If we let the leading edge of the flat plate be the origin of coordinates. 
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so X = 0 at the leading edge, and let the pressure be constant until 
Xq, where it begins to rise, the local skin friction will be zero, 
according to Stratford, in the region of increasing pressure when 

d C 

(2 Cp )^/4(n-2) -33—^ = 1.06 eClO"^ (81) 

3 is a constant equal to about .66 for Reynolds numbers on the order' 

of 10^, R is the local Reynolds number based on Uq and x, and n is 

^ 6 7 

abolit 6 or 7 for Reynolds numbers between 10 and 10 . Equation (81) 

is not sensitive to the value of n, and a change in value of n = 6 to 

n = 7 will yield about 4 percent effect. Experimental data has indicated 

that a value of n - log,« R , where R is the Reynolds number based on 

®0 ®0 

Xq and Uq, yields an error of less than 1 percent. 

Equation (81) was derived on the basis of an inner and an outer 

Tl 2 

solution. However, at the point where C = , the join of the 

inner and outer solutions reaches the outer edge of the boundary layer 
when idealized velocity profiles are used, and the equation is no longer 
valid. 


Equation (81) is a differential equation for C (x) resulting in 

^0 

the desired zero skin friction Velocity gradient, The solution of this 
equation yields 


C f—)= 0.645 {0.435 R 

Po'^V ®0 


2/n 


n - 2 


^Pn - n 1 


(82) 


Stratford extended the range of C beyond C 

% Po Po 


n - 2 
n + 1 


by assuming a 


constant shape factor, H - -g- . Using this assumption, is defined 

by ' 


C = 1 - — 

Po (x+b)^/^ 



n - 2 
n + 1 


( 83 ) 
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d C 


where a and b are constants such that C and . 

Pn dx 


Po 


are continuous at 


C ri ^ 2 

Pq = . Stratford [10] conducted a series of experiments that 

showed that equations (82) and (83) do indeed yield a near zero value 
for the skin friction at the surface. 

Rather than work with the pressure coefficient, we will work with 
the normalized velocity, q(s) . To make the conversion, we use the 
relationship 


C = 1 - fr 

Po ^ 


(84) 


and the relationship 
U 


q = 

Therefore, 

q = qr 


u "o 


U=o Uq u„ 


^0 u. 


(85) 


[i-c f^\~\ 

L Pov^o J - 


1/2 


( 86 ) 


Equation (79) can now be written as 


C = — 
^ "u 


s s 

f ' 


u 


u 


q(s)ds + 


0^ 


^0 


II - c i I 
L Po *0 J 


1/2 


ds 


(87) 


0 


We now need the relationship between x and s. Along the flat plate, 
X is the coordinate, while s is the coordinate along the surface of the 
airfoil. In order to relate the two, we set the momentum thickness 
of the boundary layer at s = s^^ equal to the momentum thickness of 
the boundary layer at x = x^ is then 
Sq'- 


"0 


I / U(s) \ j_ 

J V U, 


( 88 ) 


0 


' 0 - 
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The relation between s and x must be such that Sq and x^ occur at the 
same location. If we define the constant k such that 




the relation between x and sXbecomes 
s == X + (k-1) X 


(89) 


(90) 


or 


\ 


X = s - (k- 1) X 


0 


(91) 


We can now define a new variable z such that 


z = 


(k-1) 


0 0 

With this relation, equation (88) can be written 
1 


1 = 


q(z) I dz 


1-k 


^0 


C92) 


(93) 


We will use the notation Z to indicate the value of z at s = S: 

With this notation, equation (87) can be written as 
1 Z , 

2 


u 


Z+k-1 


{ 


q(z)dz + q^ 


k-1 


II - (z)]^^2dz} 

^0 


(94) 


The problem is to achieve a maximum from equation (94) . The 

length of the upper surface (s^) is a factor that will be fixed by the 

desired chord length of the airfoil and C is a function of R , 

% 

which is a function of the free stream Reynolds number, and is, 
therefore, a design parameter to be specified for the airfoil. 


voatsoHEiBaakwcM 
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Therefore, in order to maximize the lift, there are three variables 

in equation (94) that can be altered to maximize the lift. In the 

first integral, q(z) has no limitation other than it is non-decreasing. 

In the second term, q^, the maximum velocity on the upper surface, has 

not been specified. The third variable which is implicit in equation 

(94) is Sq, the location of the beginning of the pressure rise. 

The first integral will be maximized if q(z) is a maximum. This 

indicates that q(z) should accelerate instantaneously to the valve q^ 

and remain there until s = s^. This means that equation (93) is now 
1 

1 = I dz (95) 

1-k 

so k = 1. The instantaneous acceleration of the velocity after the 
stagnation point indicates that the stagnation point is on the leading 
edge, which results in a sharp leading edge, which will yield poor 
results in off-design angles of attack. More will be said about 
this in a later chapter. 

As can. be seen from equation (83), C reaches the value unity 
only at infinity. Since the chord length of infinity is not feasible, 
a non-zero trailing edge velocity will be accepted, which will be re- 
ferred to as q^. This leads to a sharp trailing edge, which is accept- 
able for our purposes. Therefore, there are now three values that can 
be varied to lead to a maximum G. , q„, q , and s„. These three quanti- 
ties are not independent: once two of them have been determined, the 

third is also determined once R and s are specified. Therefore, if 

®0 ^ 

the Value 6f one of these quantities is Specified, one of the other quan- 
tities may be varied to obtain a maximum G^^, and the third value will be 
uniquely defined. 
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I£ Sq is specified, for any value of a will result. However, 
as q^ increases, increases, and there is no upper limit. If qQ is 
fixed, the increases as increases, and reaches a maximum when s^ * 
s^, or there is no recovery at all, and q^ = <!,,• This is not a satis- 


u 


factory solution. The only remaining possibility is to fix q^ and vary 


q^ and s^ to get a maximum C^. 


The function C can take either the form of equation (80) or the 

Y 

form of equation (83), depending on the values of R and — . For the 

®0 ^0 

values of R normally encountered in airfoils, C will reach the value 


n - 2 
n + 1 


soon after the pressure recovery begins. The point at which this 


occurs will be referred to as s ( or z or x , depending on the coordinate 

m ^ m m’ ^ * 

system used) . 

If the notation 


a = 




(96) 


b = 


/X 


(97) 


is used, equation (83) can be written 


c CZ) "1 - -^W2 
Pq (z+b )^^ 


C < -ii— 

Po 


(98) 


At the trailing edge , z = Z 


(Z) - T - 


Pq' ' (Z+bV'^^ 


Or, using equation (86) 


*lu ' '•o 




(100) 
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This can be rearranged to 

qo - Cz + 


/~T 

a 


(101) 


Therefore, equation (94) can be rewritten as 
“■u 


2 q XZ+b')^/"^ 


m 


L Z+k-1 /IT 
2/N 


{ 1 + 


D) 


1/2 . ^ ^ 
dz + 


dz 




1 - 0.645(0.435 R 
- ®0 


(102) 


m 


Since 2 ^ depends only on /the first integral is a function only of 

R and can be abbreviated as 
®0 z 


m 


I(R^ ) = 

0 1 


1 - 0.645(0.435 R ^^^Cz^'^^-1)) 

®0 ^ 


1/2 


dz (103) 


By carrying out the integration in the last term in equation (100) and 


utilizing the fact that k = 1, 


2 q 

C, = — i- {1 t 1(R ) > i 


z rp 


®0 ^ 


cz . bV '"' 


(z -bV/‘' 

^ m 


If R and q are specified, C. varies only with z. Therefore, the 
6 q -u L 

maximum value of occurs at a value of z where the first derivative 
of with respect to z is zero. Taking the derivative and setting it 
equal to zero yields 


3/4 


1 ^ I(R ■) - I /T" (z^.b) 

C Q O III 


(z+b') + IVT' b’ (z+b') 


3/4 


[ 


1 + i(R ) - 4 (z +bV''‘^ 

^ 


=• 0 


(105) 
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’ 1/4 

Equation (105) is a fourth degree algebraic equation for (Z+b ) , 

I 1/4 

and, as such, yields four roots for (Z+b ) . For Reynolds numbers in 

the range of airfoils for General Aviation aircraft, there are two complex 

roots, one positive root, and one negative root. The only root that is 

physically meaningful is the positive real root. By looking at the 

value of the second derivativ.e, it can be shown that the positive root 

I 1/4 

does indeed yield a maximum Substituting the value of (Z+b ) 

back into equation (101) yields, q^ as a function of q^, and, by solving 

Z = rp (106) 


for Xp, Sq can be determined . Therefore, the variational problem has 
been solved. 

Up to this time, it has been assumed that the flow in the boundary 
layer has been fully turbulent from the stagnation point to the point 
of initial pressure recovery. However, in actual flows, the flow will 
be initially laminar, and will transition to a turbulent flow at some 
point beyind the stagnation point. With laminar flow present in the 
initial boundary layer, equation (88) is no longer valid, but will 


be replaced by 


V “o 

S, 


U 


"t -J 


0 


U 3 


ds 


(107) 


where the subscript t indicates the variable is evaluated at the transi- 


tion point." Equation (107) assumes an instantaneous transition with a 
continuous momentum thickness at s^. The point s,^. can occ,ur at s^, but 
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the purposes of this paper, we will assume that s^ is not greater than 
Sq, as Stratford only considered the separation of fully turbulent flows. 


If g is defined as 


(108) 


and s is replaced by z, as in the case of fully turbulent flows, equation 
(106) can be written 


1 =38.2 


^0 ^0 


+ (l-g)k 


(109) 


In the case of the initially laminar boundary layer, another value 

is required, namely R , the Reynolds number at transition. R will 

0 6 
cr cr 

be defined as 




8 =0 % 


8 ^0 % 


Note that the maximum velocity on the airfoil is used to define R , 

6 

cr 

rather than the free stream velocity. Rearranging terms, we get 


g k - 


R V R 
e e 

cr cr 

^0 ^0 ^e 


( 111 ) 


Substituting equation (110) into equation (108) yields 


k = 1 - 38.2 R R 

=0 ®cr 


For the case of flow that is fully laminar until the pressure 
recovery starts, s^= Sq and g = 1. Therefore equation (111) becomes 
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Substituting this into (111) yields 


0=1- 38.2 R R 

®0 ®cr 


(114) 


or 


R = 38.2 R 

®cr 


5/8 


(115) 


Hence, 


R 


^ R 


cr 


38.2 e. 


(116) 


This indicates that for R greater than 1.655 x 10^, R can be larger 

0 ®cr 

than R , even though s^. is greater than s . 

Once again, to find the maximum the derivative of Cj^ with 

respect to z needs to be set to zero. However, in this case, k is not 

equal to one, but is determined by equation (112). Taking the derivative 

of (102), using equation (103) to define I(R ) , and setting the result 

®0 

equal to zero yields 


3 

4 


K + I(R ) - ~ /T' Cz^+b')^/"^ (Z+bO - ? (K-l-b') (Z+b’)^/^ 

*3 m o 


j CK-l-b'). 


K + I(R^ ) - ~ (Z+b ')^/'* 


= 0 


(117) 


As in the case of fully turbulent flow, this is a fourth degree algebraic 

1/4 

equation for (Z+b') . The roots for R *s of interest here are similar 

' ' ^0 ■ 

to those of the fully turbulent case, and the single real positive root 
yields the maximum . The difference between the fully turbulent velocity 
distribution and the partially laminar velocity distribution is the location 


of the start of the pressure reco\ery. For the partially lam' nar flow, s_ 
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is greater than the Sq of the fully turbulent flow. Therefore, for 

equal , s^, and an airfoil with partially laminar flow will 

yield a slightly higher C, . 

max 
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IV. USING THE EPPLER METHOD TO DEVELOP A STRATFORD AIRFOIL 

In the design o£ an airfoil, the first parameter to he chosen is 

an approximate value for the velocity at the trailing edge. The 

higher the q^, the higher the of the airfoil, but, also, the greater 

the angle of the trailing edge. The chord Reynolds number is dependent 

on R and a , so the choice of R will be dependent on the desired 
®0 ^ 

chord Reynolds number and q^. By study of Tables 1 through 3 and with 

foresight gained through experience, an approximate R can be determined 

®0 

This value can be modified during the design process until the desired 

chord Reynolds number is obtained. 

If the airfoil is to have partially laminar flow along the rooftop, 

R will have to be chosen. One possible choice for R would be the 
e e 

cr cr 

critical Reynolds number for a flat plate, which, according to Schlicting 
. [10] is 3.2 oc 10^. From equation (112) , K can now be determined. By 


setting C = y solving equation (82) for ~, z can be 

Pq ^ ^ ^0 ^ d C 

determined, a and b are now determined such that C and Pq are 

- 2 Pq dx 

continuous at C = • Equation (117) can then be solved for 

nl/4 


n + 1 ■ 

(Z+b which determines Z. By substituting this value of Z into 


equation (101), the ratio between q^ and q^ is determined, However, the 
exact value of q^ or q^ is not kno\m until the Eppler problem is solved, 
Using equations (90) and (106), the relation between x^ and X can be 
determined as 


^0 Z + (K-1) 


where Xq is non-dimensionalizedw'ths^.. 

Since the relation between t and (p is not known e.tJ niter the 
solution of the Eppler problem, the location of the stmt of - ae pressure 
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recovery is not known in terms of the coordinates used to specify the 
input to the Eppler problem. However, using the approximation 

X = ^ + cos(j)) (119) 

an approximate value of can be obtained as 

cos({) X 2 - 1 (120) 

0 

where Sq = K Xq. 

In figure 3, if the initial slope of (w') and the value of 
at <J> = 0 are specified, the values of K and ]1 in equation (39) will 
be fixed* The value of cj)^, the start of the cusped region at the trail 
ing edge, should be chosen as small as possible while maintaining a 
reasonable trailing edge angle and airfoil thickness. An approximate 
range of values of is 2° to 30°. The upper limit corresponds to 
higher values of q^. 

The values of j w', w, and d) on the lower surface must also be 
Uif s 

specified. Since the primary interest is in the pressure distribution 
on the upper surface at the design angle of attack, these Values are 
not critical, but they can be altered to obtain better performance 
at off-design conditions. More will be said about this in the next 
chapter. In equation (39) , cannot be specified, as it will be 
determined by the solution to the problem* However, the solution 

A ^ 

can be iterated by Varying either a Or K (or K.) on the upper or low- 
er surface (or both surfaces) until = K^, a previously deter- 

mined value. This allows control over the trailing edge angle. For 
the purposes of this paper, a will be allowed to vary on the lower 
surface . 
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The next design criteria to be chosen are the angles o£ attack., All 

o£ the angles o£ attack in the Eppler problem are re£erred to the zero 

li£t line, rather than the chord line. The angle o£ attack on the upper 

surface will be the design angle o£ attack. This is the parameter that 

will control the velocity on the upper surface. The value required to 

obtain the desired approximate can be determined by study of Tables 1 

and 3, The angles of attack on the upper and lower surface must meet 

the requirements of equation (47) . 

It has been found through experience that it is impossible to 

match the initial slope of the pressure recovery of the Eppler output 

to the Stratford output with a continous angle of attack along the upper 

surface. However, by specifying a lower angle of attack in the region 

of the initial pressure recovery, the pressure distributions can be 

matched. This means the airfoil must be broken into four regions. In 

the first region, 0 = a* will be the design angle of attack. 

In the second region, = (p , oi will be less than the design 

* 

angle of attack. In the region > 'the a will again 

L 

be equal to the design angle of attack. In the last region, which 
is the lower surface, = 2tt, cx will be the angle of attack on 

the lower surface. 

Using the input data derived above, the Eppler problem is solved, 

and maximum velocity on the upper surface and s , the surface length on 

^u 

the upper surface are determined. Using these values and the R , ^/q , 
z , X„j a’, and b' determined earlier, a Stratford distribution is generated. 
By comparison the Stratford pressure distribution with the Eppler pressure 
distribution, the necessary changes to the Eppler input can be determined. 
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The primary variables to be changed are (p , which changes the location 

w 

of the pressure recovery, w, which changes the trailing edge velocity, 
w and a^, which changes tlie initial slope of the pressure recovery, 
and {j)^j which changes the pressure distribution after the initial pressure 
drop. 

These changes are input into the Eppler problem, and a new air- 
foil is generated. Since this airfoil will have a different maximum 
velocity on the upper surface and the stagnation point will be at a 
different location, changing the value of s^, a nev^^ Stratford distri- 
bution will be needed, and the cycle repeats itself until an airfoil 
is generated that matches the corresponding Stratford distribution. 
Determining how much to change each input variable to get a desired 
change in the output of the Eppler problem is an art that can be learned 
only by experience. 

Table 1 shows the variation of several airfoils with varying 

design angle of attack and varying Table 2 shows the variation 

with varying R . These airfoils are labelled by a number which gives 
®0 

some of the pertinent information about that airfoil. The first two 

digits indicate the free stream Reynolds number {x 10^) of the airfoil. 

The second pair of digits represent the 10) of the airfoil as 

indicated by the Eppler program. The third pair of digits is R /R 

,5 ®cr ^chord 

(x This is the location (in terms of s) of the transition point. 

, 1 ■ ■ 2 

The fourth pair of digits is — ^ (x 10 ), which is the free stream Mach 
number where flow on the rooftop first becomes sonic. The last pair of 
digits is the thickness of the airfoil in percent of the chord length. 

Thus, an airfoil labelled 1640-20-34-21 indicates an airfoil with a design 
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TABLE 1 

VARIATION OF DESIGN INPUTS WITH VARYING ANGLE OF ATTACK AND 


Design 


Airfoil 

angle of 
attack 

Ip 

w 

w 

°^2 


K 

s 

1640-20-34-14 

34.00 

36,73 

5.04 

10.00 

35 . 10 

45.00 

1642-20-33-15 

36.00 

37.60 

5.07 

11 .90 

35,60 

55.00 

1543-21-33-18 

36.00 

38.00 

5.13 

3.00 

36.85 

115.00 

1447-23-31-15 

40.00 

39.00 

5.14 

7.00 

37.50 

65.00 

1447-23-32-16 

40.00 

39.00 

5.18 

2.00 

37.90 

85.00 

1348-25-31-17 

40.00 

39.60 

5.19 

2.00 

38.65 

115.00 

1449-23-31-16 

42.00 

40.00 

5.15 

8.00 

38.50 

70.00 

1350-23-31-17 

42.00 

40.00 

5.18 

5.00 

38.75 

90.00 

1350-23-30-18 

42.00 

40.40 

5.22 

2.00 

39.40 

115.00 

1350-23-28-19 

42.00 

40.50 

5.26 

-0.90 

39.80 

135.00 

1351-23-30-16 

44.00 

41.00 

5.16 

9.00 

39.50 

75.00 

1352-25-30-17 

44.00 

41.00 

5.26 

6.00 

59.70 

95.00 

1352-25-29-17 

44.00 

41.00 

5.35 

3.00 

39.90 

115.00 

1253-26-29-19 

44.00 

41.75 

5.28 

7.00 

40.50 

135.00 

1254-27-29-18 

46.00 

42.00 

5.40 • 

6.00 

40.70 

115.00 

1254-29-27-19 

46.00 

42.45 

5,40 

9.00 

40.90 

135.00 

1156-29-27-18 

48.00 

43.00 

5.46 

7.00 

41.70 

125.00 

1256-29-27-17 

48.00 

43.0^ 

5,46 

6.00 

41.80 

135.00 

1256-27-28-18 

48.00 

43.50 

5.50 

9.75 

41.95 

145.00 

1157-29-27-19 

48.00 

43.70 

5.49 

9.75 

42.20 

155.00 

1157-29-28-19 

48.00 

43.80 

5.53 

9.00 

42.35 

165.00 

1157-29-27-20 

48.00 

43.90 

5,56 

8.50 

42.50 

175.00 


TABLE 2 


VARIATION OF DESIGN INPUTS WITH VARYING R 

e 


Airfoil 

R 

"o 

-e- 

w 

“2 

-e- 

' 

K 

s 

1256-27-28-18 

1.56 

43..50 

5.50 

9.75 

41.95 

145.00 

1357-25-28-19 

1.90 

43.10 

5.45 

9.75 

41.40 

145.00 

1657-20-28-19 

2.40 

42.90 

5.42 

9.75 

41.20 

145.00 

1857-18-28-19 

2.65 

42.80 

5.39 

9.75 

41.05 

145.00 

2057-16-28-19 

3.00 

42.50 

5.41 

9.75 

40.75 

145.00 

2257-15-28-19 

3.40 

42.10 

5.36 

8.75 

40.40 

145.00 

2557-12-28-20 

4.00 

42.10 

5.37 

9.75 

40.35 

145.00 

2857-12-27-20 

4.50 

42,00 

5.32 

9.75 

40.30 

145.00 

3157-10-28-20 

5.00 

41.70 

5.30 

6.00 

40.40 

145.00 

3757-08-28-20 

6.00 

41.05 

5.31 

6.00 

39.50 

145.00 

4057-08-28-21 

6.70 

41.10 

5.27 

9.75 

39.10 

145.00 

4557--07-29-20 

7.50 

41.05 

5.29 

9.75 

39.25 

145.00 

515,7-06-28-22 

9.00 

40.85 

5.22 

9.75 

38.80 

145.00 

2256-15-29-19 

3.40 

42.05 

5.33 

8.75 

40.35 

130.00 

2057-16-28-19 

3.00 

42.25 

5.41 

7.75 

40.65 

135.00 
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TABLE 3 

VARIATION OF SELECTED PARAMETERS 


Airfoil 


^0 

<^4 


k 

s 


1640-20-34-14 

1.385 

2.907 

15.31 

56.49 

3,423 

3.591 

1642-20-33-15 

1.418 

2.990 

22.30 

59.14 

4.076 


1543-21-33-18 

1.462 

3.104 

17.72 

57.73 

6.888 


1447-23-31-15 

1.493 

3.173 

26.57 

61.02 

4.835 


1447-23-32-16 

1.506 

3.219 

25.08 

60.51 

5.758 


1348-25-31-17 

1.538 

3.249 

23.43 

60.00 

7.185 


1350-23-31-17 

1.548 

3.301 

27.53 

61.57 

6.159 

4.012 

1350-23-30-18 

1.569 

3.357 

25.83 

60.99 

7.323 


1350-23-28-19 

1.588 

3.371 

24.66 

60.62 

8 . 269 


1351-23-30-16 

1.571 

3.374 

30.96 

62.97 

5.621 

4.338 

1352-25-30-17 

1.589 

3.401 

29.83 

62.57 

6.560 


1352-25-29-17 

1.611 

3.413 

28.87 

62.24 

7.505 


1253-26-29-19 

1.626 

3.470 

27.42 

61 . 75 

8.431 


1254-27-29-18 

1.650 

3»489 

31.73 

63.45 

7.694 

4.053 

1254-29-27-19 

1.666 

3.513 

30.60 

63.06 

8.626 


1156-29-27-18 

1.691 

3.580 

33.89 

64.42 

8.341 


1256-29-27-17 

1.690 

3.620 

32.94 

64.06 

8 . 782 


1256-27-28-18 

1.710 

3.598 

32.95 

64.09 

9.278 


1157-29-27-19 

1.719 

3.615 

32.37 

63.89 

9.741 


1157-29-28-19 

1,729 

3.636 

31.74 

63.66 

10.200 


1157-29-27-20 

1.737 

3.651 

31.11 

63,44 

10.670 

3.985 

1357-25-28-19 

1.696 

3.615 

32.44 

63.90 

9.242 


1657-20-28-19 

1.692 

3.605 

32.38 

63.88 

9.238 


1857-18-28-19 

1.684 

3.624 

32.04 

63.75 

9.217 


2057-16-28-19 

1.684 

3.599 

32.19 

63.81 

9.224 


2257-15-23-19 

1.673 

3.588 

31.96 

63.73 

9.203 


2557-12-28-20 

1.675 

3.571 

32.15 

63.81 

9.209 


3857-12-27-20 

1.668 

3.570 

32.01 

63.76 

9.194 


3157-10-28-20 

1.660 

3.571 

31.80 

63.69 

9.182 


3757-08-28-20 

1.648 

3.562 

31.46 

63.57 

9.156 


4057-08-28-21 

1.641 

3.569 

31.25 

63.49 

9.143 


4557-07-29-20 

1.643 

3.563 

31.33 

63.52 

9.148 

4 .140 

5157-06-28-22 

1.630 

3.571 

30.93 

63.37 

9.122 



40 


free stream Reynolds number (based on the chord length) of 1.6 million, a 
design C. of 4.0, boundary layer transition at s/chord length ■ .20, a 

It ^ 

critical Mach number of .34, and a thickness of 21 percent* , 

■ i * 

The airfoil is divided into small segments that are of equal length 
in the circle plane. For the purposes of this study, the circle was 
divided into 92 segments, but more or fewer points could have been used. 

As the number of points is increased, the accuracy of the problem in- 
creases, but the number of calculations also increases. The data relat- 
ing to a position on the airfoil (all the tj) data) is given in terms of 
these circle divisions. The relation between a vClue of (j) given in degrees 
and the same value of given in circle divisions is 


(}) Ccircle divisions) = <j> (degrees) 


number of circle divisions 


360 


( 121 ) 


The data in Tables 1 to 3 is given in this form. 

In column 2 of Table 1, the design angle of attack is listed. The 

dependence of the trailing edge velocity, listed in Table 3 column 1 

can be noted here. In column 2 of Table 2, the R of different airfoils 

®0 

is listed. The design angle of attack of the airfoils in Table 2 is 48®, 

and the R of the airfoils in Table 1 is 1 .56 million. Goljumn 3 in both 
. ®0 ' ^ ' , . . '■ 
tables lists (J)^, the start of the pressure recovery and column 4 lists 

w, the value of W at 9 =0. If the entire upper surface were at the 

' i ; 

design angle of attack, this would be the same as /q^. However, the 
discrepancy is caused by the short segment after the initial pressure 
recovery, where a* is much lower. This second angle of attack is listed 
in column 5, and is listed in column 6. The vafue of a* is «2 i*' 
the range Column 7 lists K^. which is the value 


iterates to . 
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For the remainder o£ the design inputs, the following values were 

used: 

cj) = 1.5 

w' = 50.00 

T = 30.00 

(b 

^W 

6 = 4.99 

w’ = 80.00 

R = 3.2 X 10^ 

e 

cr 

Table 3 lists some of the results of the above input data. Column 4 
lists a on the lower surface. Although this value is determined through 
iteration in the program, if the initial guess is as close as possible 
to the final result, fewer iterations are required. Column 5 lists (j)^ , 
which is referred to as the stagnation point. In actuality, this has 
no physical meaning, as it is the stagnation point when each section of 
the airfoil is at its own a^. The stagnation point of the airfoil at 
any one angle of attack will be different than (f) . 

/l .. .. 

Column 6 lists the Cj^ as obtained from the Lockheed [11] program. 
This check was not run for all the airfoils, so this list is not complete 
Of those that were run, the values of C did not agree With the C pre- 
dieted by the Eppler program. This is apparently due to the inability 
of the flow to make the sudden recovery at the trailing edge. In the 
next chapter, an improvement will be suggested that should yield better 
results. 

A typical airfoil is shown in figure 4. Figure 5 portrays the 
pressure distribution of this airfoil, and figure 6 is the equivalent 
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pressure distribution from the Lockheed program. It can bd noted that 
the rooftop pressure distribution is not nearly as negative with the 
Lockheed distribution as it is with the Eppler distribution. This is 
typical of all tbe airfoils tested with the Lockheed program. Figures 7 
and 8 show several more airfoil profiles and pressure distributions. 

In figure 4, a small protuberance can be noted at the trailing 
edge on the lower surface. This is typical of the airfoils with the 
higher trailing edge velocity, and is apparently due to the trailing 
edge angle becoming extremely large (greater than 180®} In many cases, 
this protuberance caused a failure of the Lockheed program due to the 
method of determining the chord line in this program. The Lockheed 
program starts at the trailing edge and moves along the lower surface 
until it finds a point where the distance from the trailing edge to 
to the point is less than the distance from the trailing edge to the 
previous point. This can occur in the region of the protuberance, 
v;hich leaves too few points on the lower surface for a reasonable 
solution. This problem can be eliminated by smoothing the protuber- 
ance out of the profile. By comparing results before and after 
smoothing with airfoils that did work in the Lockheed program, it 
was found that this protuberance had very little effect on the results 
of the Lockheed program* A method of eliminating this protuberance will 
be suggested in the next chapter. 
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Figure 8. Several pressure distributions 
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V. SUGGESTIONS FOR FURTHER STUDY 

There are several directions in which continued study of airfoils 
with a Stratford distribution could be channeled. Perhaps the most 
pressing is the matter of what trailing edge velocity to use. Con- 
ventional NACA series airfoils, as listed in Abbott and Von Doenhoff 
[13], have a trailing edge velocity in the neighborhood of .8Ueo, which 
is significantly less than the values attained in the present study. 
While it is generally agreed that the trailing edge velocity of an air- 
foil must be limited, a cursory search of available literature has 
not indicated what the maximum value can be. Smith [14] demonstrates 
that increasing the trailing edge velocity will increase the lift, 
but indicates that the method of increasing the velocity above .SU;,, 
is through the use of flaps, and does not discuss the possibility of 
increasing trailing edge velocity through design of the airfoil. 

Edwards [5] indicated difficulty with the thickness of the airfoil and 

I 

trailing edge shape when the trailing edge velocity exceeded 1.08. 

The trailing edge velocities of the airfoils in the present study are 
quite possibly too high, indicating that the separation of the trailing 
edge will move forward on the upper surface, destroying the pressure 
distribution. What needs to be done is, first, a systematic review of 
the literature to find if a maximum value of trailing edge velocity has 
been determined, and, second, if nothing is found in the literature, a 
maximum trailing edge velocity should be determined, either through 
experimental or analytical methods. 

Another inconsistency in applying the Stratford distribution is in 
the determining of In equation (93) (for fully turbulent flow) or 
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equation (107) (for partially laminar flow), the assumption was made 

that was equal to one from the stagnation point to the point where 
0 

the pressure rise began. In order to maximize the lift at design angle 
of attack, this is truej However, this implies that the stagnation is 
on the leading edge, which leads to a sharp leading edge. This sharp 
leading edge means a negative pressure peak at the leading edge in off- 
design conditions, which will result in separation of the boundary layer. 
What actually happens in the design of the airfoils is the stagnation 
point moves down on to the lower surface, and moves back towards the 
trailing edge some distance. For example, the stagnation point of the 
1657-20-38-19 airfoil is at about ^ = 0.18, The velocity forward of 
the stagnation point remains at a low level until after the flow has 
gone around the leading edge, where it accelerates to the rooftop velocity. 
To correct for this error, the Stratford program should be rewritten, 
using equation(93) or equation (107) to define k. 

As noted in chapter 4, a small protuberance is generated on the 
lower surface at the trailing edge, due to the large trailing edge angle. 
One method of eliminating this large trailing edge is to increase the 
length of the cusp ed region. This was done to the 1657-20-28-19 airfoil, 
with the resulting profile shown in figure 9 and the resulting pressure 
distribution shown in figure 10. The cusped region in this case was in- 
creased until the resulting x's in the airfoil profile were in monotoni- 
cally increasing order on the upper surface. The resulting was 8.0 
circle divisions, or 31.3°. It was necessary to adjust (})^ and (pj in 
order to match the Stratford distribution, but the remaining inputs are 
the same as the original airfoil. The resulting <|)^ was 42.5, compared to 
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42.9 on the original airfoil, and the resulting was 40.8, compared 
to 41.2. The geometric angle of attack changed to 22.69® from 19.06°. 

One unexpected improvement was the change in design lift coefficient 
predicted by the Lockheed program. The original airfoil had a design 
Ci of 3.868, while the airfoil with the modified trailing edge had a 
design Cj^ of 5 . 845 

The concentration of the design of airfoils in the present study 

was on the upper surface. The on the lower surface was the variable 

allowed to vary in the iteration to set Perhaps one 

area of further study could be designing the lower surface such that 

a Stratford distribution occurs on the lower surface as well as on the 

upper surface. This would result in an airfoil with better off-design 

performance, as the flow would probably remain attached for all angles 

of attack between the design angle of attack on the upper surface and 

the design angle of attack on the lower surface. If the design angle 

of attack on the lower surface was the negative of the design angle of 

attack on the upper surface, and R was the same for both upper and 

lower surfaces, the resulting airfoil would be a symmetrical airfoil. 

This design of symmetrical airfoils suggests a further possible use of 

the program, the design of low drag struts. These low drag struts would 

simply be symmetrical airfoils, with = 0.0 (the design angle of 

attack) in the region ^1— region = 

L 

would be some positive angle of attack. As is increased, the 
thickness of the strut will increase. In order to get a true Stratford 
distribution, the modification redefining k described earlier would need 
to be implemented, as it is impossible to have another constant velocity 
on the upper and lower surfaces simultaneously, 



53 


to have a constant velocity on the upper and lower surfaces simultaneously. 

Another improvement to the design method would be to combine the Eppler 
and Stratford program, so no visual comparison of the two velocity distri- 
butions would be necessary. This could be done by changing the form of 

W to be of the form 
w 



where C (-— ) is given by equations (82) and (83) . R would be an input 

PQ ^0 ®0 

parameter to the problem. Assuming ‘i’j— '^l then be equal to q^, 

and the velocity vrould be of the form of equation (86) . The form of 

could remain as it is to allow control over the trailing edge angle. It 

might then be necessary to iterate on (|)^ until equation (105) is satisfied 

in order to achieve a maximum lift. 

One of the undesirable features of the airfoils designed in the present 
study is the large positive pitching moment. This moment might be reduced 
with a sacrifice of some lift by allowing the suction rise at the leading 
edge to occur more slowly. This would result in a less negative pressure 
region in the forward region of the upper surface, and thus a smaller mo- 
ment. This slower suction rise could be input by specifying to be 
greater than the design angle of attack in the region (p . If k 

were modified as described earlier, the value of k would be increased 
by this slower suction rise (since ^ would be less than one over much 
of the rooftop region) , causing the start of the pressure rise region 
to move farther aft, regaining some of the lift lost by the loss of 
the low pressure region at the leading edge. 
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The airfoils of the present study were all designed with a critical 
Reynolds number of 3.2 x 10^. However, data [15] from tests of the 
University of Illinois HL^ 1720-00 airfoil indicates the flow remains 
laminar throughout the rooftop region. The results of the Lockheed [121 
program tend to support this result on the airfoils designed in the 
present study. Therefore, perhaps for future studies the critical 
Reynolds number should be defined by equation (116). 
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APPENDIX A 
THE EPPLER PROGRAM 

The calculations required for the solution of the Eppler problem 
are carried out with the aid of an IBM-360/75 computer at the University 
of Illinois. The Eppler program, not only determines the profile of the 
airfoil but also determines the boundary layer momentum thickness and 
the energy form parameter. However, in the present application, the 
boundary layer capabilities of the program have not been fully utilized. 

The required programs are kept in files on the PLORTS system. 

The file name of the Eppler program is EPPLER, while the file names 
of the required input data are EDATA through EDATG. A sample imput 
data deck is shown in Figure kl, 

The first card in an input data deck for the Eppler program is 
a card with an Alpha-numeric listing of the titles of the cards that 
follow. .These titles are read in 2 0A4 format. It is essential that 
the order of the titles not be changed and all titles must be included 
on this card, even if the named card is not used in the program. This 
first card can be thought of as part of the program itself, as it is 
never changed. The remaining cards, with the exception of the title, 
are in the format (A4,16,14F. 2) , Some of the data that is input through 
the P5.2 format is divided by a factor of 10 in the program, so it is 
important not to specify the decimal point. All the data should be 
right justified, and the program will convert the data to the correct 
multiple of 10. The data that is divided by 10. in the program will be 
identified in the following discussion as having a psuedo-format of F5.3 
The manner in which the data on each card is treated is determined 
by the title, which is listed in the first 4 spaces on each card. The 



TRA1TRA2ALFAAGAMABSZ REENDEBETAPLOTTITL 
ABSZ 9200 100 

AGAM 100 100 100 100 100 100 

TITL 

U OF I HLE 1657-20-28-19 AIRFOIL 

PLOT 3470-83372128815625330711089424000 

TRA10000027 4120 4800 4290 0975 0000 4800 9200 3200 

TRA2000027 150 4290 100 5000 0542 3004000 10080000499 20014500001 

BETA -1 100 

RE 03 01633 

ENDE 
/* 


Figure Al. A sample input data deck 
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data is read into the program as MARKE, NUPU, and PUFF, where MARKE 
is the title, NUPU is an integer, and PUFF is a 14 element array. The 
data is then transferred to the appropriate variable according to the 
title. 

The first title listed on the first card is the TRAl title. The 
TRAl card is the card that inputs the (|)^ and a^. The (J)^ are input in 
terms of circle divisions, and the (j). are input in degrees, is 

determined by the program, so it is input as zero. The (f)^ and are 
input as pairs, and up to seven pairs can be input on one card. If 
it is desired to break the circle plane into more than seven segments, 
more TRAl cards need be specified; with a maximum of four cards, as 
storage is allowed for only 28 segments. The last must be equal 
to the number of divisions in the circle. The d). must be listed in 
increasing order, including the computed value of . 

Spaces 5 through 10 of the TRAl card (NUPU) are reserved for 
the profile number. If several different airfoils are developed at 
the same time, they can be identified by this profile number. 

Spaces 5 through 10 of the TRA2 card are also reserved for the 
profile number, but in this case, the profile number is used only to 
keep track of the input data, as this number is not used in the program 
These spaces can also be left blank on the TRA2 card. 

The remainder of the words on the TRA2 card define the input 
Velocity function. Words 1 through 5 define the upper surface and 
words 6 through 10 define the lower surface. Word 1 is (j)^, given in 
circle divisions, and word 2 is The meaning of words 4 and 5 

depends on the word 3. If word 3 is 0.0, word 4 is k and word 5 is p. 
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\ 

I 

I£ word 3 is 1.0, word 4 is w and word 5 is w. If word 3 is 2.0, 
word 4 is lii and word 5 is w. Words 4 and 5 are divided by 10.0 in 
the program, so the psuedo format is F5.3, 

The specification of w and w' Cword 3 being equal to 1.0) is 

• . 

recommended only with large values of w , so the path of W is strongly 

w 

curved. The process converges slowly when w' is small, and convergence 

is not guaranteed when ]i is negative. For less strongly curved paths, 

the specification of p and w is recommended Cword 3 equals 2.0). 

Words 6 through 10 define the lower surface in the same manner 

that words 1 through 5 define the upper surface. Thus, for a symmetrical 

airfoil, words 6 through 10 would repeat words 1 through 5. 

Word 11 is referred to as ITMOD, and determines the variable that 

is changed in the iteration process to set to the specified value. 

The ITMOD is 0.0, no iteration is carried out. If ITMOD is 1.0, the 

on the upper surface are altered by a factor Aa^ until attains 

the desired value. If ITMOD is 2.0, the a. the lower surface will be 

1 

altered and if ITMOD is 3.0, the will be altered on both the upper 
and lower surface by an equal amount. If ITMOD is 4.0, K is modified, 
if ITMOD is 5,0, K is modified, and if ITMOD is 6.0, K and K are modified 
by equal amounts. ITMOD = 3.0 or 6.0 is useful for symmetrical airfoils. 

Word 12 is K^, written in the psuedo -format of F5.3. Word 13 is 
the tolerance acceptable in the computation, also written in the pseudo- 
format of F5.3. A suggested value for this is ,001, the smallest value 
available in the F5.3 format. Word 14 is not used. 

The next card in the list is the ALFA card. This card inputs the 
various angles of attack that the pressure distribution is developed for 
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and that are used in the boundary layer portion of the program. The first 
word after the title is NAL, the number of angles of attack listed, in 16 
format. NAL can be as large as 14. If NAL is specified as larger than 14, 
it is reset to 4. The next 14 (or less) words are the angles of attack, 
in degrees, written in F5.2 format. If NAL is given as a negative number, 
the angle of attack will be given on the TRAl card, where i is on the 
ALFA card in F5.2 format (see the sample data deck in Figure A1 for an 
example of this) If an ALFA card is given with no angles of attack and 
NAL=0, the angles of attack of the previous profile are repeated. 

The AGAM card controls the output of the Eppler program. The 16 of 
the AGAM card is ignored, but 14 AGAM(i)'s are read in F5.2 format. In 
general, the AGAM(i) 's are either zero and not zero. If AGAM(l) is not 
zero, the X and y coordinates of the airfoil are generated. If AGAM(l) 
is equal to zero, only the transcendental equation is solved. If AGAM (2) 
is not equal to zero, the profile list will be printed, along with a 
velocity distribution for each angle of attack on the ALFA card. If 
AGAM(3) is not zero, the input data and the solution to the transcendental 
equation is printed out for the initial input and the final iteration. If 
AGAM(4) is not zero, the input data and the solution to the transcendental 
equation will be printed out for all iterations. AGAM(5) and AGAM(6) refer 
to the boundary layer portion of the program. If AGAM(5) is not zero, the 
program will print out a listing of the distance along the surface from 
the stagnation point, the local velocity, the energy thickness form para- 
meter (the energy dissipation boundary layer thickness divided by the 
momentum thickness), and the momentum thickness. If AGAM (5) is equal to 
1.0, the local Reynolds number, based on the momentum thickness and the 
local velocity is printed out instead of the momentum thickness. If AGAM(6) 
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is not equal to zero, the boundary layer transition point, boundary 
layer separation point, and drag (calculated by the Squire Young Method) 

are printed out. AGAM(7) through AGAM(14) are not presently used, but 
are reserved for further use. 

At the University of Illinois, most runs are made with AGAM(l) 
through AGAM (6) equal to 1.0. This results in the most complete output. 

An attempt to run with AGAM(6) equal to zero resulted in the failure of 
the program for unknown reasons . 

Card ABSZ lists the number of circle divisions, NKR, in spaces 11 
through 15. NKR must be divisible by 4, and NKR + 1 points result in the 
profile of the airfoil. As NKR is increased, the accuracy of the solution 
increases, as well as the computational time required. The maximum NKR is 
120, but 60 is usually a sufficient number unless large slopes in the velocity 
function are encountered, as with a Stratford distribution. For the airfoils 
designed at the University of Illinois, an NKR of 92 was chosen. 

The ABSZ card also lists ABFA in spaces 16 through 20, which multi- 
plies all values given in circle divisions. ABFA is normally equal to 
1.0. It is necessary to change ABFA only if the number of circle divisions 

is changed, so it is not necessary to change all the input data given in 

circle divisions. If no ABSZ card is given, NKR is set to 60 and ABFA is 

set to 1.0. 

The RE card is used to input the Reynolds number into the program. 

The pseudo-format of the RE card is (A4, 6X, 5(211, 3X, F5,5)). The first 
of the II words represents MA, which at one time was used to determine the 
suction mode. Since the capability of boundary layer suction has been re- 
moved from the program, this word is no longer used. The second II word 
is MU, the mode for boundary layer transition. When MU is equal to 1, 
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transition is by laminar separation. If MU is equal to 2, transition occurs 


at the first decrease in velocity. If MU is equal to 3, transition occurs 


when the velocity remains constant throughout a step distance or decreases. 


If MU is 4, transition occurs vdien the natural logarithm of the local Reynolds 


number based on ^2 local velocity exceeds or equals 18.43 H ^2 - 21.74. 


MU = 5 is similar to MU = 4, except the value that In (RE) is compared to is 


18.43 H 22 - 22.10. Therefore, MU = 5 is a more conservative estimate for 


transition. The F5.3 word is the free stream Reynolds number, based on the 


chord length and free stream velocity. All lengths in the program are non- 


dimensionalized with respect to this chord length, and all velocities are 


non-dimensionalized with respect to this velocity. There can be up to 5 


Reynolds numbers, each with its own MA and MU. The program will continue 


to read in Reynolds numbers (up to 5) until a zero value is read as a 


Reynolds number. 


The ENDE card is necessary for proper termination of the program. 


It is the final data card, and indicates all data has been read in, 


The next three titles on the list are cards that have been added to 


the program at the University of Illinois. The first of these cards is the 


BETA card, which replaces the ALFA card. If a BETA card is used instead of 


an ALFA, card, either a punched output is generated or data is filed into 


the PLORTS system that is used by the Stratford program. This data con- 


sists of four parts, written in 6F 12, 9 format. The first part is DS, the 


increment of the surface distance for each x increment. There are NKR OS's 


generated. The other three parts are a velocity function (VF) , and x and 


y coordinates of the airfoil. There are NKR + 1 of each of these values. 


The velocity function is equal to the local velocity divided by (1 + Cosa^) . 
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The program was originally designed to give a punched output, but was 
modified to file the data directly into PLORTS. However, as the PLORTS 
system is due to be removed from the IBM-360 at the University of Illinois, 
it will be necessary to change back to a punched output deck. 

The next card that has been added to the program is the PLOT card. 
This card reads data into the system that is then either punched out or 
filed into PLORTS. Nothing is done with this data by the program, as this 
is only a convenient method of getting data into the input deck for the 
Stratford programs. 

The last card to be described is the TITL card. No data is on 
the TITL card, but this card signals that the next card is in 20A4 format, 
and is the title of the airfoil. This title will be printed in the output 
and inserted into the Stratford input deck. 

There are some restrictions on the order the cards are read in. 

The ABSZ (if one is used) , AGAM, TITL, and PLOT cards should be read into 
the computer first, although not necessarily in that order. The data 
on these cards remains valid until another similar card is read into the 
computer. Thus, for example, if several profiles are to be developed 
with the same number of circle divisions, it is not necessary to repeat 
the ABSZ card. The next cards to be read in are the TRAl and TRA2 cards, 
in that order. Once the TRA2 card is read in, the profile is generated. 

The ALFA or BETA card is then read in, followed by the RE card. The RE 
card initiates the calculation of the boundary layer. If other profiles 
are desired, new TRAl and TRA2 cards can now be read in, preceded by new 
ABSZ, AGAM, PLOT, and TITL cards, as necessary. These cards can be followed 
by ALFA or BETA and RE cards if boundary layer inform^-tion is desired. 
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The ENDE card terminates the program after all the profiles and boundary 
layer calculations are complete. 

The descriptions of the output which follows assiames AGAM(l) through 
AGAMC6) are not equal to zero. If any of these words are equal to zero, 
the corresponding portion of the output will be deleted. 

The first data listed in the output are the input data and the 
solution to the transcendental equation. This data is preceded by the 
title, profile number, iteration number, and iteration mode (0 through 6). 

The headings of the table of data do not agree with the nomenclature 
presented in this paper. NUE represents the same quantity as i|)^, ALPHA 
is a^, WS is w and w, iVHK is w and w', DRAK is K and K, DRAM is y and 
y, HK is and K^, FLA is (j)^^ and and LAS is <})^ and cj)^. 

The next data listed are the profile of the airfoil in x and y 
coordinates and the velocity distribution for each angle of attack on the 
ALFA or BETA card. AT the end of this listing, the values of CM, BETA, 

ETA, SX, and SY are printed out. CM is the moment coefficient at zero 
lift and BETA is the angle between the zero lift line and the chord line. 
Since all angles of attack are given in reference to the zero lift line, 
this angle is necessary to compute the geometric angle of attack. ETA, 

SX, and SY are apparently remnants of trouble shooting the program, as 
they are not particularly useful. ETA is the number of points in the circle 
plane divided by the chord and Tt. This term is used in non-dimensionalizing 
the chord. SX and SY are summation of the x and y coordinates of the airfoil 
profile. 

The last section of data is derived from the boundary layer portion 
of the program. First there are two tables, one for the upper surface and 



64 


one for the lower surface. These tables list the surface coordinate , local 
velocity, '^2’ AGAM(S) is equal to 1.0, the local Reynolds 

number based on 6 ^ and the local velocity is printed in place of ^ 2 . 
However, nothing in the output indicates that this has been done, so it • 
is important that it be noted that AGAM(5} is equal to 1.0 if this data 
is to be used. If “32 ^ negative number, the flow in the boundary 

layer is turbulent, 

Follcwxng these two tables are listings for the upper and lower 
surface transition points, separation points, and drag coefficients. 

Once again, there is a problem of nomenclature, as the transition points 
are under the heading INS., the separation points are under the heading 
TRANS., and the drag coefficient are under the heading SEP.. The transi- 
tion and separation points are given in terms of surface coordinates. 

The plotting routine for the Eppler airfoils is filed in two 
PLORTS files, PLOTMN and PLOTOBJ. The data- for the plotting routine 
is normally filed in PLDATCN) and PLOTOBJ (N), where (N) is a number 
between 1 and 5, The first card in PLDATCN) is a card of the form 
b^NAMlbN=92,ALPHA=18.42,$END where b is a blank space, N is the number 
of points on the airfoil (normally 92) , and ALPHA is the design angle 
of attack. This card is a punched output card of the Eppler program, 
but it is not the first card. Therefore, the deck must be rearranged 
to be in the proper order. The next part of the PLDATCN) file consists 
of the DS, VF, X, and y cards, as punched out by the Eppler program. 
Finally, the first card is repeated four times, but with different 
angles of attack. The only data that is changed in the PLOTOBJ (N) 
file is the first card, which is the title card in 2 0A4 format. The 
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remainder of the data in the PLOTOBJ(N) file is concerned with numbering 
the axiSj and always remains the same. 

The plotting files are run in the following order: PLOTMN, PLDAT(N) , 
PLOTOBJ, PLOTOBJ(N). The first (or main) portion of the program (that 
part filed under the PLORTS file PLOTMN) determines the pressure coefficients 
circulation, lift coefficient, and center of circulation, first for the 
design angle of attack and then for the other four angles of attack listed 
on the last four cards in FLOAT (N) . The airfoil coordinates and pressure 
coefficients are stored on tape. The second part of the program (PLOTOBJ) 
then runs, reads the data on tape, and plots the airfoil profile and 
pressure distribution for the design angle of attack. 
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APPENDIX B 

THE STRATFORD PROGRAM 

The Stratford program is divided into two parts. The first part, 
which is kept in the PLORTS file TABLE, determines the roots of equation 
(117). The second part, kept in the PLORTS file STRPLOT, takes the out- 
put data from TABLE and the Eppler program and plots a Stratford distri- 
bution that corresponds to the RED input to TABLE with an initial velocity 
equal to the rooftop velocity of the Eppler airfoil. 

There is no external input data for the program in the file TABLE. 
If a different set of data is desired, the changes have to be made in the 
program itself. Therefore, for example, the statement RECR - N.NEN, 
where N.NEN is the desired critical Reynolds number, must aopear early 
in the program. The program is set up to solve equation (1 l 7) for up 

to 30 values of R . If less than 30 values are desired, the statement 
®0 

NUM=30 must be altered to reflect this. The first value of R is input 

®0 

through the statement RE0(1)-N.NEN, where N.NEN is the desired value. 

The remaining values of R are input through the statement RE0(K+1)= 

REO (K) +N . NEN , where N.NEN is the desired step size. The trailing edge 
velocity can be input through the statement VTE=NN.N, but, since all the 
data except Xq and the chord length are non-dimensionalized, this value is 
of no consequence. If a value of kinematic viscosity (v) other than 
1,6 X 10~^ is desired, the statement ANU=160.H-6 can be changed to 
reflect this . 

The roots of equation (117) are determined through the use of a 
subroutine from the IBM Scientific Subroutine Package [17J named POLRT. 
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This subroutine determines all of the roots of the polynomial, and the 

program searches through these roots until it finds the positive real 

root. From this real root, the value of Z can be determined. Through 

the use of equation (101), the ratio of can then be determined, 

and, from equation (106) , the value of can be determined. 

For each value of R , the TABLE program prints out values of k, a', 

®0 

b', D = ^u/qj^, Z, z^, Xq, and the chord length. The values of a', b’,/ 

D, z , Z, k, and R are read into the Eppler program through the PLOT 
m ®Q 

card, and are then output with the rest of the Eppler output either on 

cards or filed into PLORTS. This data is used by the second part of the 

Stratford program, which is filed in the PLORTS file STRPLOT. 

The program in STRPLOT takes the data from the Eppler program, 

and, through the use of the Calcomp plotter, draws the required Stratford 

pressure distribution based on the data from TABLE and the Eppler program. 

The first card read into STRPLOT is the title, written in 20A4 format. 

The second card contains N, SFl, and SF2, in 14, 2F10.7 format. N is a 

number that is no longer used in the Stratford program, and can be left 

blank. SFl is a scale factor in the x direction. In order to match the 

output of the Eppler plot program^ this should be 10.0. SF2 is the scale 

factor in the y direction, and should be 0.31'25 to match the Eppler plot. 

The next data read in are a', b', D, U„, z , Z, s , and k, in 8F10.7 

U m u 

R is then read in F 15. 4 format. The last two sets of data read in are 

®o ; 

the X coordinates of the airfoil and DS, the surface distance between 
the points on the airfoil. Both of these sets of data are in 6F12.9 format 
The program is presently set up for 93 points on the airfoil . 
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APPENDIX C 

THE LOCKHEED PROGRAM 

The Lockheed program was used as a method of checking the results of 
the Eppler program. Given the coordinates of an airfoil, the Lockheed 
program determines the lift, drag, and moment coefficients. The theory 
and application of the Lockheed program is documented in references 12 and 
17. However, there have been a few modifications to the program as run at 
the University of Illinois. These modifications will be the subject of 
this appendix. 

The input cards to the program are the same as in reference 12 ex- 
cept for cards 2 and 3. Card 2, which is concerned with the plot subroutine 
that is not used, was eliminated. Card 3 has two more variables, IPLOT and 
MXTRAP. IPLOT is presently not used, but is reserved for use in conjunction 
with a plotting subroutine. MXTRAP will be explained in the following pages. 

The major modification to the program was the restoration of the 
smoothing process of the local Mach number at the trailing edge. As noted 
in reference 18, large Mach number gradients at the trailing edge create 
undesired "kinks" in the equivalent airfoil camber line. In order to correct 
for this, the computed Mach numbers at the last two points bn the upper sur- 
face of the airfoil are discarded, and a linear least-squares fit is applied 
to the last five remaining points on the upper surface , The least squares 
fit is then shifted until a smooth fransition occurs at the most forward 
point (i.e,, at the seventh point from the trailing edge) , and the curve is 
extrapolated to replace the last two points. The computed Mach number at 
the last three points on the lov^er surface is then discarded, and a second 
order interpolation between the last poinf on the upper surface and the 
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fourth and fifth points from the trailing edge on the lower surface is 
used to redefine the last three points on the lower surface. These 
modified values of local Mach niunbers are used only to determine the 
local boundary layer characteristics. The actual computed values are 
printed out in the output of the Lockheed program. 

When North Carolina State University modified the multi-element pro- 
gram to obtain the single element program described in reference 12, they 
found that the smoothing and extrapolation resulted in a significant lift 
from symmetrical airfoils at zero angle of attack. Therefore, they removed 
this portion of the program. Studies at the University of Illinois have 
shown that correlation with analytical results (at least for airfoils 
with, a Stratford distribution) is better with the smoothing and extrapola- 
tion routine in the program. Therefore, the smoothing and extrapolation 
routine was restored, but vath two modifications. First, a second order 
least squares fit was used on the upper surface, and, second, the number 
of points smoothed on the upper surface was made a variable called MXTRAP. 

If MXTRAP is 0, no smoothing is done, while if MXTRAP is a positive integer 
greater than 2, this number Of points are smoothed on the upper surface. 
MXTRAP should be at least 3 so a second order least square fit can be 
done. Most of the work at the University of Illinois has been done with 
MXTRAP=7, so the smoothing routine is the same as the original multi-element 
Lockheed program . 

For each angle of attack, the resulting lift coefficient is divided 
by the drag coefficient . This lift over drag data is then printed out in 
the table at the end of the output, along with the lift/ drag, and moment 
coefficients . 
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The punch option was modified to make it compatible with the plotting 
routine at the University of Illinois. The first punched output card is 
the title, in 20A4 format. The next card contains the reference chord 
length, stagnation temperature, chord Reynolds number, Prandtl number, 
heat transfer coefficient, MXTRAP, and the number 1 in (5F12. 5, 215) format. 
This last number 1 indicates to the plotting routine that the single element 
program was used, as opposed to the multi-element program. The third card 
contains the number of elements in the airfoil (always one) and the number 
of points in the airfoil (always 65) in 215 format. The fourth card lists 
the number of free stream Mach numbers and the number of angles of attack 
in 215 format. The next sets of data are the x and y coordinates of the 
airfoil in 6F12.8 format. 

The remaining data is repeated for each angle of attack and each 
Mach number. The first two cards list the separation points for the upper 
and lower surface. Often, due to an unknown problem in the program, an 
extra card is punched out at this point, indicating no transition on either 
the upper or lower surface. Therefore, before the output deck can be used, 
it must be checked to make sure there are only two cards listing separation 
points for each angle of attack and Mach number combination. The next 
card lists the free stream Mach number and angle of attack in 2F12.5 format. 
This is followed by the lift coefficient, drag coefficient, and moment 
coefficient in 3F12, 5 format. The last portion of data is the local pressure 
coefficient at each of the 65 points on the airfoil in 6F 12. 8 format. 

This punched output is. then fed into the plotting program, which is 
filed in PLORTS under the file name PICT. This program- plots (using the 
calcomp plotter) , first, an outline of the airfoil profile, with a listing 
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of the, reference chord length, stagnation temperature, chord Reynolds nvimber 
Prandtl number, and heat transfer coefficient. The program then plots a , 
p'ressure distribution for each angle of attack and Mach number. If at 
least three angles of attack have been specified, the program then plots 

■ I ’ ; : ; 

a C, versus angle of attack curve, a versus C. curve, and a C versus C. 
curve. If the boundary layer separates at some point before c = 0.95, 
the point of separation is indicated on the versus angle of attack curve. 
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